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Regular two-component bouncing cosmologies and perturbations therein 
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We present a full investigation of scalar perturbations in a rather generic model for a regular 
bouncing universe, where the bounce is triggered by an effective perfect fluid with negative energy 
density. Long before and after the bounce the universe is dominated by a source with positive 
energy density, which may be a perfect fluid, a scalar field, or any other source with an intrinsic 
isocurvature perturbation. Within this framework, we present an analytical method to accurately 
estimate the spectrum of large-scale scalar perturbations until their reentry, long after the bounce. 
We also propose a simple way to identify non-singular gauge-invariant variables through the bounce 
and present the results of extensive numerical tests in several possible realizations of the scenario. 
In no case do we find that the spectrum of the pre-bounce growing mode of the Bardeen potential 
can be transferred to a post-bounce constant mode. 



I. INTRODUCTION 

Motivations for studying bouncing cosmologies come 
from today's widespread belief that a consistent theory of 
quantum gravity should cure the initial cosmological sin- 
gularity -as well as other singularities that plague Gen- 
eral Relativity (GR)- by a suitable modification of its 
short-distance behaviour. In the context of string the- 
ory, bouncing cosmologies of different kinds emerge as the 
most reasonable alternative to the Big Bang singularity. 
In the pre-bounce phase, as one goes to past-infinity, the 
universe is driven towards a non-singular attractor where 
generic initial conditions can be chosen. 

The first proposal of a bouncing string cosmology was 
the so-called Pre-Big Bang scenario 0, which exploits 
the non-minimal coupling of the dilaton in superstring 
theory to drive a Pre-Big Bang super-inflationary epoch, 
ending when higher-order derivatives and/or loop correc- 
tions become non-negligible. A second realization of the 
same idea was given by the ekpyrotic/cyclic scenario 0, 
inspired by Horawa-Witten braneworlds 1^. Here the 
pre-bounce is characterized by the slow approach of two 
parallel branes, which eventually collide and then move 
away, giving rise to an ordinary expanding universe on 
each of the two branes. Both scenarios, when viewed 
in a four-dimensional Einstein frame, can be represented 
by a universe bouncing from contraction to a standard 
expansion. 

A crucial test for any cosmological model lies on its ca- 
pability of accounting for the origin, the growth and the 
spectrum of the fluctuations needed to explain the large- 
scale structures and the anisotropics in the cosmic mi- 
crowave background. It is well known that observations 
require a nearly scale-invariant spectrum for the Bardeen 
potential. Such a spectrum can be naturally provided by 
the initial quantum fluctuations of the inflaton in the 
standard inflationary picture Q . Alternatively, it can be 
provided by the initial isocurvature perturbations of an- 



other field, as in the "curvaton" mechanism j^, which is 
naturally implemented through the axion field in the Pre- 
Big Bang scenario 0. A third possibility, the claim that 
the ekpyrotic scenario could directlty provide (without 
appealing to the curvaton mechanism) a scale-invariant 
spectrum of adiabatic density perturbations, has given 
rise to a considerable debate about the correct way to 
describe the effects of a bounce on cosmological pertur- 
bations 7] . 

The main point at issue here concerns the Bardeen 
potential: it certainly grows very large, with a red or 
scale-invariant spectrum, during the pre-bounce phase. 
Does this mode survive the bounce or does it decay 
away, allowing other modes to take over? The use of 
classical matching conditions for perturbation through a 
space-like hypersurface (which assume the continuity 
of the comoving curvature Q is not trivial in the case of a 
bounce, since the Hubble rate changes sign and the null 
energy condition (NEC) is necessarily violated for some 
time during the bounce. Then, more exotic matching 
conditions cannot be excluded, depending on the physics 
of the bounce 0]. 

In this situation, the only possibility to test some gen- 
eral arguments is to construct explicit models of regular 
bounces where the final spectrum of the Bardeen poten- 
tial can be explicitly evaluated. In the context of GR, the 
more conservative possibility is to study closed universes, 
where the bounce is possible without violating the NEC 
^d\. In this case, the presence of the curvature scale is 
responsible for an explicit momentum dependence in the 
transfer matrix from pre- to post-bounce perturbations. 

On the other hand, if one wants to stick to spatially 
flat cosmologies, a very good approximation presumably 
after a long pre-bounce epoch, one has to introduce a 
negative energy-density component that determines the 
transition from a contracting universe to an expanding 
one, dominated by a standard component, such as ra- 
diation. In this respect, Peter & Pinto-Neto lljj stud- 
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ied a bounce induced by a ghost scalar field, which con- 
nects a radiation-dominated contraction to a radiation- 
dominated expansion. Surprisingly they found that the 
Bardeen potential keeps the spectrum of its pre-bounce 
growin g m ode. Similar conclusions have been found by 
Finelli in a model with two perfect fluids (one with 
negative energy density). 

On the other hand, explicit bounce examples have been 
built using only one scalar field and triggering the bounce 
by higher derivatives and loop corrections to the Einstein 
equations 13J or by potentials involving integrals on spa- 
tially flat sections of a compact universe [3| . Both these 
one-component examples pointed at a complete post- 
bounce decay of the Bardeen pre-bounce growing mode. 

Coming back to two-component models, Allen & 
Wands [T^ found again a decay of the Bardeen pre- 
bounce growing mode in a cosmology dominated by a 
scalar field with a tracking potential (which allows it to 
mimic the background evolution of any kind of fluids), 
with the bounce triggered by a ghost scalar field. The 
reasons why the conclusions of this work are in contrast 
with those of Ref. .1 1 remained unclear. 

In this paper we would like to bring more clarity in 
the investigation of regular bouncing cosmologies by in- 
troducing a rather generic cosmological model, which in- 
cludes, as particular cases in a much wider parameter 
space, all previous models containing two non-interacting 
components. Wc find that the prc-bouncc growing mode 
of the Bardeen potential always matches a decaying mode 
in the post-bounce, leaving no trace in the final spectrum 
at horizon re-entry. Our analysis proves that the claims 
in Refs. |]3 are incorrect, while it confirms the re- 
sults of Ref. [13 , generalizing to a much larger class of 
two-component bouncing cosmologies. A short summary 
of our results was anticipated in Ref. . 

In this work we provide all the details of the analy- 
sis summarized in Ref. |16) . In particular, we explic- 
itly provide the steps necessary to produce the analytical 
estimates from the method sketched in Ref. IQl. More- 
over, we present the numerical studies of the cosmological 
perturbations in the whole parameter space, testing and 
supporting the analytical estimates. Finally, we expand 
the discussion on the gauge-choice issue and the way to 
construct non-singular gauge-invariant variables. 

In Sect. 2, we introduce the class of backgrounds we 
want to consider, write down all the equations, and dis- 
cuss the general behaviour of our background quantities. 
In Sect. 3, we introduce scalar perturbations, discuss 
the gauge issue, introduce the most convenient gauge- 
invariant quantities for our analysis, and discuss the vac- 
uum normalization of perturbations that provides the ini- 
tial conditions for our analysis. In Sect. 4, we describe 
the method used for analytical estimates and provide all 
the necessary steps for a systematic approach of the prob- 
lem. In Sect. 5, we discuss these estimates throughout 
the parameter space that describes our bouncing cos- 
mologies. In Sect. 6, we numerically check all the es- 
timates discussed in Sect. 5. Finally, Sect. 7 contains a 



discussion of the results. 



II. A GENERAL MODEL FOR 
TWO-COMPONENT BOUNCING 
COSMOLOGIES 

A common criticism against bouncing cosmologies trig- 
gered by a ghost scalar field -or negative energy matter- 
is that they deal with unphysical sources that are plagued 
by quantum instabilities and/or inconsistencies. How- 
ever, we must keep in mind that the "true" bounce should 
be the outcome of a short-distance modification of GR, 
such as those coming from String/M-theory. Once we 
choose to perform an effective description in the Einstein 
frame, we cannot avoid (in a spatially flat cosmology) vi- 
olations of the NEC in the sense that the short-distance 
corrections must enter as an effective negative energy (or 
pressure) contributions to the right-hand side of Fried- 
man's equations. Unless this is accepted, one must give 
up using the tools of GR and wait for the full String/M- 
theory to tell us how scalar perturbations behave during 
the high-energy phase. 

In this paper we assume that a General Relativity de- 
scription is indeed possible (at least for modes that are 
far outside the horizon during the bounce) and that the 
high-energy corrections can be modelled as an effective 
new negative-energy density field pf,, which becomes neg- 
ligible as soon as we go far from the bounce itself. At this 
point it makes no sense to worry about the quantum in- 
stabilities of such an effective field, since the quantum 
description of high-energy corrections can only be done 
in the full mother theory. Nevertheless, nothing prevents 
the investigation of classical perturbations with the stan- 
dard tools of GR. 

The class of models we investigate in this paper have 
the following characteristics: 

• The universe is filled by two components a and b, 
with energy density pa > and pf, < 0, respectively. 

• Far from the bounce, pa dominates the evolution 
of the universe, while pb is important only at the 
bounce and becomes subdominant elsewhere. 

• The two components do not decay and interact only 
gravitationally. 

• Far from the bounce, the ratios 

rm=K„/Pm: m^a,b, (1) 

are asymptotically constant, that is we have a sta- 
ble attractor in the past (implying also w,„ — 

Pm / Pm — Tm)- 

The Friedman equations for a two-component cosmol- 
ogy can be simply written as 

= a^p (2) 
n'^ + 2H' = -c?v (3) 
p'-|-3H(p + p)=0, (4) 
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where a prime denotes the derivative w.r.t. conformal 
time 77, a is the scale factor, Ti. — a' /a is the conformal 
Hubble rate, and we are using units such that SttG = 1. 
The total energy density is the sum of the individual 
energy densities of the two cosmological components p = 
Pa+Pb, and the same happens for the pressure p ~ Pa+Pb- 
Since the two fluids have neither interactions nor decay, 
they satisfy the continuity equation Q independently. 

As a consequence of the constancy of r,„, the attractor 
in the past is described by the following power laws 



q 



n 



Pa 

Pb 

q = 



\v 



|2+2g 

3g2 



|,^|3g(l+r6 

2 



i + 3r„ 



(5) 
(6) 

(7) 

(8) 

(9) 



Here we have already chosen a convenient normalization 
for the scale factor and for the conformal time 7], so that 
the two fluids have comparable energy densities at 77 ~ 
— 1. Of course, the condition pb/pa — > for rj —>■ —00 is 
satisfied if and only if > Fq. 

The power-law solution holds until ?; ^ —1. Then the 
bounce occurs at 77 = 0, when Ti. changes sign and finally 
at some time 77 ~ 1 the post-bounce expanding power-law 
solution is established with the same Ta and Ti, as in the 
pre-bounce, so that Eqs. (0)-® hold again. In practice, 
we allow the bounce to be asymmetric within the range 
1 7^1 < 1, but require the asymptotic power laws to be the 
same before and after the bounce. Actually, this final 
assumption could be easily relaxed, but we do not know 
any explicit example of two-components regular bounce 
with different pre- and post-bounce power laws. In Fig. ^ 
we show typical behaviours of the background functions. 



III. SCALAR PERTURBATIONS 

The perturbed line element, including only scalar per- 
turbations, is 

ds^ = a^{Tj){{l + 2(j))dTf ~ 2B^,dr]dx' 

- [(1 - 2^)6^j + 2E^,j] dx'dx^ } , (10) 

while the perturbed energy-momentum tensor of each 
source reads: 



6pr, 



-{Pm +Pa)Vm,j 



™ " \{P7n+Pm)V7n,i - {Prn + Spm)Sij + 

(11) 

with TO — a,b respectively; V,n is the velocity potential 
and Srn is the anisotropic stress. It is worth recalling 
the effects of a gauge transformation on each of these 
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FIG. 1: Typical behaviours of the background functions in 
a bouncing cosmology. At the top we have the scale factor, 
in the middle the conformal Hubble parameter H, and at 
the bottom the function (7i^ — Ti.') = a^{p + p)/2, which 
becomes negative when NEC violation occurs. The borders 
of the bounce phase have been put in evidence by vertical 
dashed lines at the times 77 = ±1. 
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variables. For an infinitesimal coordinate change 



(12) 



the scalar perturbations are modified in the following 
way: 



4> = 




(13) 


^ = 




(14) 


B 




(15) 


E = 




(16) 


5Pi = 


- sp^ - p[e 


(17) 






(18) 






(19) 



A. A general model for source perturbations 

In order to close the system of equations, for each com- 
ponent we need to specify a relation between the pertur- 
bative variables. Since we want to make a single treat- 
ment for both perfect fiuids and scalar fields, we need a 
general enough relation, containing both kinds of sources 
as particular cases. We therefore impose, for the domi- 
nant source, 

dpa=cl5pa + aipa+Pa)Va, a = , (20) 

Pa +Pa 

while for the secondary component we choose perfect 
fiuid relations 



5pb = Cf,(5p&, = Tb- 



(21) 



We also assume the anisotropic stress to vanish for 
both components. 

In general, we allow Ta and to be different, opening 
the possibility of a non-vanishing a. In order to better 
understand the physical meaning of and a we can go 
to the comoving gauge with respect to the dominant fluid 
(Va = 0). We then have 



5pa\v^=0 
5pa\Va=o' 



(22) 



which identifies Ca as the sound speed of the a-fluid. In 
fact, this quantity is properly defined in the frame where 
the fluid is at rest. We will see later that, indeed, is 
the velocity of propagation of fluctuations in the a-fiuid, 
thus explaining our somewhat unconventional notation. 

The perfect fluid relation is recovered when the equal- 
ity Ta — Ca holds. In fact, in this case, we simply have 
a — and the pressure perturbation is just proportional 
to the density perturbation in any gauge with the same 
coefficient Ta = c^. On the other hand, if Ta ^ c^, we 
have a 0, implying the presence of an intrinsic isocur- 
vature mode 

TaSSa = Spa ~ {p'a/p'a)SPa = a{Pa + Pa) {Va ~ Spa/p'a) ■ 

(23) 



In particular, a scalar field with a self-interaction 
potential V{ip) has 

P^-^+Vi^)^ = (24) 

Spa - - 0^ + V'{^)S^, Va - 4' (25) 



6pa = ^ V - 0^ - V\ip)dip, 



(26) 



while the covariant conservation of the energy- 
momentum tensor gives 



ip" + 2nip' + a^V'{ip) = 0. 



(27) 



It is immediate to prove that Eq. (|20l) holds with = 1 
and a = —2a^Vip/p'. Notice that a free scalar field is 
completely equivalent to a perfect fiuid with stiff equa- 
tion of state, p = p, while the presence of any potential 
inevitably introduces an intrinsic isocurvature mode. 

From this discussion, we conclude that Eq. (|20|l is 
general enough to contain perfect fiuids, scalar fields and 
any source with an intrinsic isocurvature mode generated 
according to Eq. (|23|l . This ensures our treatment of a 
wide generality and, in particular, it allows us to include 
all previously studied models as particular cases. 



B. The gauge issue 

Sooner or later, all authors studying the problem of 
the evolution of cosmological perturbations through the 
bounce have met a non-trivial difficulty in making a con- 
venient choice of the gauge. The commonly used lon- 
gitudinal gauge makes the Bardeen potential appear di- 
rectly as one entry of the perturbed metric. Typically, 
this variable becomes very large in the approach to the 
bounce, making the validity of the perturbative approach 
questionable. Other commonly used gauges introduce ex- 
plicit divergences at the NEC violation or at the bounce. 
Under some circumstances, a safe gauge is provided by 
setting to zero the perturbation of a scalar field whose 
time derivative never changes sign [l4l ITsI l . 

Our approach is to assume the existence of a regular 
gauge where all the components of the metric and of the 
energy-momentum tensors stay finite. In order for the 
linear theory to be valid, the perturbations should also 
stay small. If no regular gauge exists, then we have to 
conclude that the background is not a stable solution 
of Einstein equations and thus it is uninteresting as a 
cosmological model. Our approach will be applicable to 
all regular bounces that are stable solutions of the Ein- 
stein equations. In practical terms, a regular gauge must 
have (/), B^i, E^ij, dpm, Spm and {pm + Pm)Vm finite 
and small (m = a,b). Notice that it is the combination 
{Pm +Pm)Vm that enters the energy-momentum tensor of 
the TO-source. This means that Vm can possibly diverge 
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when {pm + Pm) goes to zero provided the above prod- 
uct stays finite. The same happens for the total velocity 
potential: 



(28) 



which may diverge at the NEC violation, the product 
(p + staying finite. 

Using the regular gauge, it is easy to construct gauge- 
invariant combinations that are regular and can thus be 
safely used to follow perturbations through the bounce. 
Indeed, the Bardeen potential 



^lf = ^ + H{E' - B) 



(29) 



is a regular combination of regular perturbations. Old 
matching conditions suggest to use the comoving curva- 
ture perturbation C, = tp + "HV. However, the presence of 
the velocity potential does not protect this variable from 
possible divergences. We propose to use instead 



(30) 



which is guaranteed to be regular, since 2(7i^ — 7i')V = 
a^{p + p)V is finite in the regular gauge. We then intro- 
duce individual gauge-invariant variables 

Cm = ^a^(Pm +Pm)(V' + UV,n) (31) 



n 



which are also regular by construction. Of course 

c = Ca + d. 



(32) 



(33) 



Before proceeding with the perturbation equations, 
it is worth commenting briefly on the commonly used 
gauges. Using the transformations H13|I ~ H19|I . we can pass 
from the (assumed) regular gauge to any other gauge. 
These transformations may either explicitly induce diver- 
gences, or ensure that the variables in the transformed 
gauge are regular as well. Let us analyse some cases ex- 
plicitly. 

The zero-curvature gauge {^p — 0), obtained by the 
regular gauge setting ^ = —tpH'^, introduces terms with 
TC^^ in 0, B and Vi, which diverge at the bounce. 

The comoving gauge (V = 0), obtained by setting — 
V, introduces V in all other variables. As said before, V 
can diverge at the NEC violation. 

The synchronous gauge {(p = B ~ 0) is obtained by set- 
ting — a^^ J acpdr] and £, — J Bdrj + J drja^^ J acpdr]' . 
This transformation induces no divergences and so we 
can say that if a regular gauge exists where all pertur- 
bations are finite, then the variables in the synchronous 
gauge are also finite. 

Similarly, for the longitudinal gauge (i? = i? = 0), 
obtained by setting £, — E and = i? — no divergence 
is introduced. 



The uniform density gauge {Sp = 0) is obtained by 
setting = 5p/p'. By the continuity equation, we see 
that the transformed variables have to diverge both at 
the NEC violation and at the bounce TL — Q. 

Concerning a possible "uniform pressure gauge" {5p — 
0, obtained by setting — 6p/p') it is difficult to say 
something general. However, it is hard to imagine a 
model where p' does not vanish at any time during the 
bounce. 

The uniform individual density or pressure gauges have 
a similar fate. They typically introduce divergences at 
the bounce. The gauge comoving with one of the com- 
ponents, instead, introduces divergences only if Vm itself 
diverges, but this is not necessarily true. The regularity 
of the gauges in Refs '13, IT3. ITEI | can be understood in 
the light of this remark. 

Among those just described, the synchronous and lon- 
gitudinal gauges are the only ones that keep all the vari- 
ables finite. Yet we know that in the longitudinal gauge 
i/) = 5" and the Bardeen potential grows very large in 
the approach to the bounce. Indeed, the finiteness of all 
variables does not assure their smallness. In this case, 
one can still consider the variables in the longitudinal 
gauge as useful auxiliary mathematical quantities: the 
transformation to longitudinal-gauge variables has to be 
regarded as a simple mathematical change of variables, 
without the physical meaning of a coordinate transfor- 
mation. Nevertheless, nothing prevents one from solv- 
ing the equations in the longitudinal-gauge variables and 
then discuss the results with gauge-invariant quantities, 
whose physical interpretation is clear. In fact, the lon- 
gitudinal variables have the important property of being 
regular (even if large) and thus can be used to go through 
the bounce safely. 



C. Perturbation equations 

Let us start by perturbing the Einstein equations. The 
00, iO, M, i ^ j components respectively read 

a^5p + m{n(t) + i}') + 2P[ii! + n{E' -b)] = q (34) 

2{n(t> + ip') = a^{p + p)V (35) 
-]-a^5p + T-L^ct) + 2H> -f -I- 2^^-' + ip" (36) 



- + iUiE' -B) + {E" - B') = 0. 



(37) 



We can use Eq. to ehminate 0, Eqs. if^ and lf?T)) 
for 5pm, Eqs. and to ehminate 5pm and Vm, 

and finally Eq. H29|l for ip. Equations 00, ii and i ^ j 
become 

{iC + kH)n + ^a + '^b = <^ (38) 
C' - (a - n)~Qa - WCfc + cl^a + cl^b (39) 



, 2H-n' c 

Xj/' -I \^ = — 



(40) 



where we have made use of the fact that Ca + Cb — C- 
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The perturbations of the covariant conservation of the 
energy-momentum tensor of the first fluid, after using the 
Hamiltonian constraint (jSHlj give 



_ a^{Pa +Pa) z 



2H 



(41) 



, ^cjn^ + a^{p+p + pa+ Pa) 2 Pa + Pa 



Ca = -3a 



2n 

2 Pa + Pa;: 



a -(42) 



The perturbations of the covariant conservation of the 
energy-momentum tensor of the second fluid give the 
same equations with a ^ b. The sum of Eq. (|41|l with 
its b counterpart yields Eq. (|39|l . while the sum of Eq. 
(|42|1 with its counterpart yields a linear combination of 
Eqs. gOI) and 

Summing up, Eqs. (|41|l and (|42|) and their b counter- 
parts are an independent set of coupled first-order dif- 
ferential equations, which completely describe cosmolog- 
ical perturbations in our problem in terms of the regular 
gauge- invariant variables (^a, Cbi and ^!b- The fact 
that Ti. appears in some denominators does not indicate 
that the differential equations break down at the bounce. 
Actually, one can easily check by the definitions H31|l and 
(|32|l that the numerators also contain an Ti, factor multi- 
plied by finite functions. 

Later on, we shall widely use another set of variables 
to describe the system, namely (, Cb, ^' and '^b- This set 
has the advantage of including directly 5" and C, which 
are the most interesting quantities characterizing cosmo- 
logical perturbations. In terms of these variables the set 
of differential equations becomes 



^, 2H^-H' C 
f'H * - — 

n H 



n 

aHPb + Pb) j^2,j. 



Cb + 



2n'-n'~ a^pb+Pb);: 



H 



2H 



(43) 
(44) 

(45) 
(46) 



bitrary choice. In fact, we can think that the first com- 
ponent is an ordinary perfect fiuid, a scalar field, or any 
other field initially in a vacuum state. On the contrary, 
on the basis of the former discussion about the nature of 
the secondary component, wc cannot really justify giv- 
ing a well-defined quantum initial state to this effective 
component. Yet, in order to compare our results with 
other works, which have systematically chosen the vac- 
uum normalization even for the secondary component, 
we will follow this attitude as well. In the conclusions 
we will comment about the generality of our results and 
their (in)dependence on the initial state of the secondary- 
component perturbations. 

Taking the derivatives of Eq. H41|) and eliminating '^a 
and ^E"^ using Eqs. and we get a second-order 
equation for Ca coupled to Cb and ^'b- Of course, we also 
get an analogous second-order equation for Cb- Introduc- 
ing the Sasaki-Mukhanov canonical variables 



v„ 



= Cm/[CmiPni +P„i)^/^7Y], 



(47) 



the two equations can be written in a very compact way 
as 

+ {elk' ~ zllZa) Va = Oipb/Pa)'/' (48) 
< + (elk' - 41 Zb) Vb = 0{pb/pa)'^^ (49) 

with Za = a^Pa+Pa^^'^/iCan) and Zb = a(i+3c=)(^^ ^ 

PbY^'^/cb- Inside 0{pb/ PaY^'^ we include all the terms 
where Va, Vb appear multiplied by background functions 
going to zero asymptotically faster than \ri\~'' and terms 
where v'^, v'b are multiplied by functions going to zero 
faster than All these terms become irrelevant in 

the asymptotic past and do not affect the initial vacuum 
normalization. We also see that multiplies the Lapla- 
cian of Vm and thus represents the velocity of propaga- 
tion of the fluctuations in each component, conflrming 
our statement in Sect. 2. 

The asymptotic solutions of Eqs. H48|l and (|49|l in the 
power law regime are 



Va= VWlCaHilHcaklvl) 

«b= VW\CbH^^Hcbk\ri\), 



(50) 
(51) 



where Hj)'^ is the Hankel function of the flrst kind, Ca 
and Cb are pure numbers of order 1, and 



D. Vacuum normalization 



(52) 
(53) 



Choosing initial conditions for the perturbations of the 
two components is by no means trivial. While for the 
dominant component we can choose ordinary vacuum 
normalization, for the secondary one this is quite an ar- 



Now that we have established the initial conditions for 
the Mukhanov variables, we can derive those for Ca and 
Cf, from Eq. (gTj). For "if a and "^b we can use Eq. (jlT|) 
and the corresponding equation for Cj- 
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IV. EVOLUTION OF PERTURBATIONS 

In this section we shall estimate the behaviour of per- 
turbations before and after the bounce, giving a complete 
discussion of the final momentum and time dependences 
in terms of the parameters of the cosmological model. 

The first step is to rewrite Eqs. 143|) - (|46ll in their 
integral form 



H 



ci(fc) 



C = a^{Pa+Pa) 



C2{k) 



Cdv 



(54) 



cin 



Cbdr] ■ 



a^{Pa +Pa) 
a^{Pa +Pa) 



(55) 



n 



cm 



+ 



Tti 



n 



-Cbdr] 



Ch^a'+^'^nipb+Pb) 



Ciik) 



a^+^-ln{pb+Pb) 



2^2 



(56) 

^bdrj 
(57) 



In these equations, each variable is expressed as the 
sum of the solution of the associated homogeneous dif- 
ferential equation, carrying an arbitrary constant Ci(fc), 
plus several integrals of the other variables, representing 
the coupling terms. This structure is particularly useful 
to study perturbations outside the horizon as the ones 
we are interested in. Indeed we shall limit out attention 
to modes with fc ^ 1, i.e. to modes which exit the hori- 
zon long before the bounce and reenter much later. For 
these modes, horizon-exit can still be described by the 
asymptotic decoupled solutions l|Sn|l and (|?T|l . 

As a result, we can match the expansions (I54ll - (|57() 
to the long- wavelength limit of the asymptotic solutions. 
As a practical example, let us examine C, which far from 
the bounce is dominated by Ca- Through Eq. H5U|) and 
(I47|l . the vacuum normalization provides the following 
expression in the long- wavelength limit for this variable 



1^1 



-l-2g 



(58) 



where all the background functions have been replaced 
by their power-law behaviours ©-(O. 

Equation (|55|l in the power-law limit gives 

C ~ [c2{k) + Ci{k)\v\''^' + co(fc)] , (59) 

where we have used Eq. (|54|l to evaluate the integral of 
^' and we have kept only the lowest-order terms in krj. 
The integration constant co(fc) can be simply absorbed 
in the definition of C2{k). 



At this point, in order to match the two expressions 
(|58|l and (|59|) . it is sufficient to impose 



ci(fc) C2ik)^k-'\ 



(60) 



where we have discarded any numerical factor of order 1, 
as we shall do throughout our analysis. In the same way, 
matching C.b^ we obtain 



C3(fc) - k-"", Ci{k) ~ k"". 



(61) 



Defining F(ry) = , C, ^I^f,, C,bj , we can write the sys- 
tem ()54 |I - H57|I in the compact form 

F(77) = (??) + B{r,) J Aii])F{r,)dr], (62) 

where 

^(0)(^) =(5^,«^(Pa+Pa)c2(fc), 

-^C3{k),a^+^'^"H{pb+Pb)c4{k)^ (63) 

is the vector of the homogeneous solutions, B{ri) and 
A{ri) are two matrices that can be easily read from Eqs. 
(I54|l - H57(l . All these quantities are just combinations of 
the background functions. Then we can construct the 
full solution of Eq. (|62|l recursively 

oo 

F(77)=F(°)(7y)+^FW(7y) (64) 

1=1 

F('){r]) ^ B{t]) [ A{r])F^'-^\T])dr]. (65) 



The whole problem is thus reduced to the calculation of 
a sufhcient number of terms in the sum (|64|1 . which are 
reduced to multiple integrals of background functions. 
In the following subsections we shall illustrate how to 
calculate these integrals in the pre-bounce, during the 
Bounce and after the Bounce. 



A. Pre-bounce evolution 

As long as the power-law behaviours ©-linjl remain 
valid, the evaluation of all integrals can be easily done. 
In fact, a generic integrand f{r]) can be expanded as 
Ci{k)\r]\'' giving 



Ci{k)\ri\''dr] Ci{k)\r] 



s+l 



rj < -1. 



(66) 



Calculating iteratively the integrals, we can safely stop 
at F^^) (r]) since the next terms would just give higher and 
higher powers of fc|77| multiplying already existing terms 
(recall that we are interested in the evolution of modes 
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outside the horizon). Keeping now the lowest order in ki] 
of each mode, we find 

k-"" |,^|~9(i+3=') + Cak'"' Ivl^^-^" (68) 

Cb - fc"" |,y|l-39(l+-?) + fc--» |,^|-9(l+3c^) 

fc-^6|^|--<?(l+3c^) ^_;,^.|^|-l-2g^ (70) 

In this expression and in the following subsections, we 
introduce the factors 



Ca = 
Xa = 



k^r]'^ a — 
1 a ^ 

A;2 a = 

1 a ^ 



(71) 
(72) 



in order to write a single expression for the solutions that 
holds both when a ^ and when a = 0. 

In Eq. 1)67(1 , we explicitly see the presence of a growing 
mode in 5* carrying a red spectrum that becomes scale- 
invariant in the limit of slow contraction q ^ 0. These 
pre-bounce solutions are valid until the energy density of 
the secondary fluid becomes comparable with the dom- 
inant one. At that point the power-law solution breaks 
down and the bounce phase starts. As we have fixed this 
time to = —1, by evaluating Eqs. l|(77|l lf7?H) at ry ~ — 1, 
we obtain the order of magnitude of the different modes 
at the onset of the bounce. Since the only scale left in 
this evaluation is k, we deduce that the mode with the 
reddest spectrum always dominates at the bounce. This 
does not mean that the same mode will also dominate 
in the post-bounce, since, as we shall see later, many of 
these modes decay after the bounce. 



B. The Bounce 

We are not very interested in the specific evolution of 
all variables during the bounce. It is just interesting to 
observe that all homogeneous solutions ^(^^(r;) can be 
easily followed through the bounce as they are directly 
expressed by combinations of background functions. 

One thing that may worry the reader is the fact that 
some integrals contain diverging functions of the back- 
ground. Actually the in the integrals of ^' and C,b 
in Eq. 1(56(1 and in the integral of ^'h in Eq. 1(57(1 are 
immediately compensated by Ti. factors in ^E*, Ch and 5'^ 
inside the integrals. The solution is a bit more subtle 
for the integrals of ( in Eqs. ((54(1 and ((57(1 . For exam- 
ple, the integral of C in Eq. 1(54(1 contains a 7i~^, which 
diverges as |?7|~^. The integral thus diverges as 



but this divergence is compensated by the Ti factor out- 
side the brackets. The result is that the limit of 5* as 
0^ is finite. The same happens for the integral of 
C in (|S7)l. Then Eqs. (|31jl - lf37l) confirm the regularity of 
all variables through the bounce. 

As we shall see later, we are interested in the integrals 
covering the whole bounce phase 



fiv)d.r]. 



(73) 



The /c-dependence factors out of the integral, and we are 
also integrating out the 77 dependence. We deduce that, 
in the absence of other scales that can enter the result, 
this integral must return a value with the same order 
of magnitude of /(—I), which has been evaluated for all 
functions in the previous subsection using the pre-bounce 
power-law limit. 

For the two integrals that diverge at the bounce, we 
can replace the integral l(73() by its principal value 



lim 



f{ri)dri + / fir])dT] 



(74) 



This expression is consistent with the continuity of the 
gauge-invariant variables. As the two divergences com- 
pensate, the dimensional argument can be safely applied 
as before. 



C. Post-bounce evolution 

In the post-bounce, the evaluation of the homogeneous 
solutions i?*-*^) {r]) is elementary, since they are directly 
expressed in terms of the background functions. In par- 
ticular, we may notice that the pre-bounce growing mode 
of with the red spectrum ci(fc)fc^^ comes from the ho- 
mogeneous solution, which becomes a decaying mode in 
the post-bounce. This does not automatically mean that 
this mode becomes irrelevant in the post-bounce. We 
need to show that no growing mode develops the same 
spectrum and that there is always at least one mode that 
dominates on this decaying mode at horizon re-entry. 

In order to evaluate the integrals in the post-bounce, 
we split their integration domain into three pieces, re- 
spectively covering the pre-bounce, the Bounce and the 
post-bounce phase 



f{v)dv = / f{v)dv + / .f{v)dv + / f{v)dv- (75) 
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These integrals give 

-1 1 

J Ci{k)\'q\''dj-i Ci{k), j f {ri)d-q ^ Ci{k) , 

-1 

■n 

J cdkM'drj^cdk)[v'^' -I]. (76) 
1 

In fact, the first integral is Eq. (|66|l evaluated at 77 = —1. 
The second integral is Eq. (fz!^ . which gives a result of 
the same order of magnitude of the pre-bounce contri- 
bution. The last integral involves background functions 
in the post-bounce and can be safely evaluated using the 
power-law limit ij^)-© again. Summing up, each in- 
tegral gives a constant contribution plus a growing or 
decaying one, depending on the value of s. If s < — 1 
the integral is dominated by the constant contribution, 
which comes from regions close to the bounce and from 
the bounce itself, while if s > —1 the integral is domi- 
nated by the growing contribution at late times, far from 
the bounce. Note that any constant contribution, once 
generated, sums to the original homogeneous mode and 
follows the same evolution. 



D. Post-bounce solution in a simplified example 

In the former subsections we have explained how to 
calculate each single integral once the integrand has been 
specified. To construct the full solution, we need to pro- 
ceed step by step, evaluating the integrals recursively. To 
explain this procedure, it is better to start with a simpli- 
fied system. 

Let us consider the subsystem (, without couplings 
to the secondary component. The starting point for the 
construction of the full post-bounce solution is given by 
the homogeneous solutions. We thus set 

*W = ^ ci(fc)fc-2^-i-2« (77) 

C^°^ = aHpa + Pa)c2{k) ^ C2(fc)77~'- (78) 

As said before, the homogeneous modes survive un- 
changed during the bounce and after the bounce. 

The first step is to evaluate the integrals using the 
homogeneous modes 

J a'^{pa+Pa) 

cmr^ W'^" - A ■ (80) 

Every integral gives a time-dependent contribution and 
a constant one (taking also into account pre-bounce and 
Bounce contributions, which are of the same order). 



As a second step, we calculate the integrals using ^I^^^' 
and C^i) 

J a^KPa+Pa) 

C2(fc)fcVM^' -!]• (82) 

In the same way we can build '^'^^^ and q'^\ and so on. 
The complete solution is the sum of all the partial con- 
tributions 

oo 

* = (83) 

oo 

C = E^~''^- (84) 

However, since we are interested in the solution just 
before horizon re-entry, we can drop all terms that are 
just k'^rf' times some other term already present in the se- 
ries. In fact, these terms determine the first oscillation of 
the mode at horizon re-entry. For example, the first term 
in the brackets of 4"^^^ and C^^^ are just k'^rf times '1''°' 
and i^^^-' respectively. The last terms in the same brack- 
ets are just fc^ times '^'^^^ and C'"' respectively, so that 
they are just k^ corrections to the homogeneous modes 
(remember that we are studying modes with k -^V). In 
this simplified example, all ^^'^ and C*^*^ with i > 3 give 
subleading terms according to this prescription. We are 
finally left with the following post-bounce solution 

fc-s+'Jyy-i-a'J + fci-? (85) 

where we have replaced Ci{k) by their explicit spectral 
dependence. 

For each variable, the first term is the homogeneous 
mode, the second and the third come from the inte- 
gral of the homogeneous mode of the other variable, the 
fourth term is the only new term coming from the sec- 
ond recursion. Now we can discuss the evolution of these 
four modes in the post-bounce. We can distinguish two 
cases; fast contraction-expansion {q > 1/2) and slow 
contraction-expansion [q < 1/2). 

Let us start by the case q > 1/2. Just after the bounce 
(?7 ^ 1), the dominant mode is always the reddest one. 
We see that is by far dominated by its homogeneous 
mode, while ( is dominated by the integral of 4". How- 
ever, when we evaluate the modes at horizon re-entry 
(77 ^ 1/fc), C is still dominated by the third mode, which 
comes from the integral of the homogeneous mode of ^' 
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and in particular is the contribution coming from regions 
close to the bounce. On the contrary, at 77 ~ \/k the ho- 
mogeneous mode of \1/ has decayed below the fourth mode 
in Eq. (|85|l . This constant mode comes from the integral 
of the dominant mode of This mode has the original 
spectrum of the homogeneous mode of multiplied by 
fc^, which is typically blue. 

In the case of slow evolution {q < 1/2), both 5* and 
C are dominated by their homogeneous modes at the 
bounce (ry ~ 1); C does not acquire a mode dominating its 
homogeneous mode and, consequently, the fourth mode 
of ^ remains subdominant till re-entry. In this case, the 
dominant mode in ^ at re-entry is still the homogeneous 
decaying mode, endowed with its red spectrum, which be- 
comes scale-invariant in the limit g ^ 0. Of course, this 
situation cannot be compatible with observations, since 
the potentially scale-invariant spectrum is carried by a 
decaying mode that becomes subdominant with respect 
to a constant mode just at horizon re-entry. Moreover, 
as we shall see later, the interplay with the secondary 
component (which we have neglected in this simplified 
treatment) is essential in this regime, since it provides 
those isocurvature modes that dominate at horizon re- 
entry. 

The two cases of this simplified example are synthe- 
sized in Fig. 12 where we have plotted 5* and C, (which is 
more familiar than C to people working in cosmological 
perturbations). Following [13, we have chosen Log 0(77) 
as time coordinate. We see that in Fig. where we have 
chosen a model with q > 1/2, ^ grows in the pre-bounce, 
dominated by the integral of Then it sets to a con- 
stant value in the post-bounce. The Bardeen potential 
pre-bounce growing mode turns into a decaying mode in 
the post-bounce, which becomes subdominant when ^' 
becomes of the order of In Fig. I^b, where we have 
chosen a model with q < 1/2, ( ~ 77^ C is dominated by 
its homogeneous (constant) mode in the pre-bounce and 
in the post-bounce. The decaying mode of ^' becomes 
of the same order as the constant mode just at horizon 
re-entry. In this case, for very small q, we have an al- 
most scale-invariant spectrum for the decaying mode of 
^, while its constant mode has a blue spectrum. 

Notice also the presence of two spikes in C before and 
after the bounce. These are the mentioned divergences 
at the NEC violation times, which make ( unsuitable for 
the description of the evolution of scalar perturbations 
close to the bounce. 



V. PERTURBATIONS AT THE ONSET OF 
HORIZON RE-ENTRY 



(a) 




Log a 



' — ' ' — " — ^ — ^ — I ' — ^ — ^ — ^ — ' ' — " — ^ — Loq a 
-1.5 -1 -0.5 0.5 1 1.5 2 



FIG. 2: (a) Simplified perturbations in a bouncing cosmology 
with a = 0, cl = 0.4, cl = 2.5. The solid line is Log]*! and 
the dashed line is LogjCl. (b) The same as (a) with — 1.6. 

all gauge-invariant variables. 

After four recursions in the integrations of Eqs. H54|l - 
(|57|l . no new modes are added to the post-bounce solu- 
tion. Its final form reads 

Cak^'ll^^^" ^Cak''^ (87) 

C~ k^'^ri"'^ + k'''^7^-^-'^i + 

k-'-r]-^ + Cak-^-r^-^-^'i + 

Xak-"-?]-^ + Cak-'"'?]-^-^'^ + 

Cak'^^r^ + Cak'^^r^-^^ (88) 



In this section, we apply the method explained in the 
previous section to the whole system of integral equa- 
tions governing our problem. We shall thus build the 
full solution of scalar perturbations in the post-bounce, 
keeping all the relevant modes. We shall then find out 
the dominant modes at the onset of horizon re-entry for 



fc-<>i,^-l-9(l+3c?) _^ ^-^,+2^-2g ^ 
^.,+2 l-g(l+3c?) ^ fc'>.+2 -2, 



(89) 
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(90) 



For each mode, we have just written their k depen- 
dence, given by the constants Ci{k), and their 77 depen- 
dence, discarding any numerical factor. As for the pre- 
bounce solution, we have used the factors Ca and Xa to 
include the case a = 0. Let us explain how these factors 
arise. The coefficients Ca come because of the integration 
of C,b in Eq. H55|l . If a vanishes, this term is absent and 
the same modes come only through integration of 5'^, 
which contain, however, an extra k'^rf' factor. Finally, 
the Xa coefficients come from subtle cancellations in the 
solution for two perfect fluids aa = 0, which require some 
more discussion. 

The first cancellation is in the integral of in Eq. 
(|55|l . In fact, the integrand is an odd function of the back- 
ground. In this situation, the pre-bounce, the Bounce 
and the post-bounce contributions cancel exactly, leaving 
no trace except for a decaying mode in the post-bounce. 
Successive recursions do not involve odd functions. The 
terms reappear multiplied by a fc^, but with the same 
time dependence. This determines the Xa factors in Eqs. 
(113 and PI). 

The second cancellation involves several integrals and 
can only be demonstrated using the exact solution of the 
two-fluid background (available in terms of hypergeomet- 
ric functions). One has to pick out the integrals of in 
^'fc and C (the latter corrected by the contribution of 
to C) and plug them into Eq. (|57|l . where the cancellation 
occurs. This explains the Xa in Eqs. H89|) and (|90|l . 

Comparing Eqs. l|S7J|-linni) to Eqs. ^-^IQ, we see 
that each mode in the pre-bounce gives rise to two modes 
in the post-bounce, one with the time-dependence of the 
homogeneous mode and one with the time dependence 
of the integral of the partner homogeneous mode. In 
some cases, these modes may be multiplied by k'^rf' or 
k"^ factors. 

As in the simplified example of the previous section, 
also here ^ has no constant mode with a possibly scale- 
invariant spectrum. At most, if there are no other modes 
dominating, the scale-invariant spectrum is carried by a 
decaying mode. This is hard to reconcile with experi- 
mental data. Nevertheless, we shall prove that there is 
always another dominant mode. 



A. Discussion of the case q 7^ 

All modes essentially depend on the characteristics of 
the background, which are encoded in the 3 parameters 
4, Fa (or equivalently g = 2/(1 -I- SFq)), eg = F;,. How- 
ever, comes always as a numerical coefficient and never 
enters the exponents of k and ?/. Following our policy of 




FIG. 3: In the plane {Va,ci) we show the interesting region for 
bouncing cosmologies in the case a 7^ 0. The shaded regions 
are excluded as giving no initial contraction, no bounce, or 
bad vacuum normalization of the secondary component. 



discarding all numerical factors, we can discuss our re- 
sults in terms of the remaining two parameters. 



Of course, not all values of F^ and eg are relevant to 
bouncing cosmologies. In fact, in order to have a bounce, 
the ratio Pbj Pa must grow during contraction. This re- 
quires eg > F^. Moreover, F^ should be greater than 
— 1/3, in order to keep g > and have an ordinary con- 
traction. Finally, eg must be greater than zero for a cor- 
rect vacuum normalization. The allowed region is shown 
in Fig. El 



In principle, we should also exclude all values of eg > 1, 
since they cause fluctuations in the secondary component 
to propagate faster than light. However, since we are 
mainly interested in the effects of the kinematics of the 
bounce on the perturbations, in order to explore the slow 
evolution region q < 1/2, we shall also include the region 
eg > 1 in our analysis. 



Comparing all the modes in Eqs. H87|I - H9U() . we can dis- 
tinguish four regions depending on the mode dominating 
each variable. These regions, shown in Fig. ^ and their 
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0.5 1 1.5 2 2. 



FIG. 4: In the case a 7^ 0, we have four regions. In regions 
yli, Ai the spectrum of the dominant component prevails in 
'5, while in regions B\, B2, the spectrum of the secondary 
component prevails. 



dominant modes are 



B2 



^ r. 

*6 



-2 

-l-2g 



- ^-''''+277-29^ 



77 



77 

77-1-2?. 



(92) 
(93) 
(94) 



The borders between the different regions are 





-A2 : 


4 


_ 1 + 

2 


(95) 


A2 


-Si : 




= 1 


(96) 


Bi 


: 




= 1. 


(97) 



We see that ^ is dominated by the spectrum of the 
primary component for all bouncing cosmologies where 
the secondary component has c2 < 1 (regions Ai and 
A2). Perturbations are also adiabatic, since C ^ 772c is 
constant outside the horizon. On the contrary, for c2 > 1, 
including also slowly evolving cosmologies (Fa > 1), is 
dominated by the spectrum of the secondary component. 
For Fa > 1 (region ^2)7 this mode is isocurvature, while 
for Fq < 1 (region Si), this mode is adiabatic. The 
distinction between regions Ai and A2 comes from the 
different time dependence in and ('b- 

Allen & Wands studied a bounce where the principal 
component is a scalar field with a tracking potential, so 




FIG. 5: In the case q = 0, we have four regions. In regions 
A and AB the spectrum of the dominant component prevails 
in ^ , while in regions B\, B2, the spectrum of the secondary 
component prevails. 



that c~ = 1, but Qfa 7^ and Fq can assume any value 
between —1/3 and 0. Their secondary component is a 
free ghost scalar field (c2 = 1, a;, = 0). Their case thus 
lies just on the border line between region A2 and Bi. 
From Eqs. H52|) and H53|l . i^a = I'b and the spectra of the 
two components coincide. They find that the final form 
of ^P, after the homogeneous mode has decayed, is the 
one given in Eqs. H92|l or (|93|l . 



Discussion of the case q = 



The case F^ = and F^ = cl contains all bounces 
developed by two perfect fluids or free scalar fields. In 
this case it is possible to write analytic solutions of the 
background in terms of hypergeometric functions, which 
are very useful to check for particular cancellations of 
modes, as explained before. In fact, all bounces described 
in this class enjoy time-reversal symmetry and thus the 
integrals of odd functions vanish. Therefore, the classi- 
fication presented in the previous subsection has to be 
revised in order to take into account these cancellations. 
Furthermore, we cannot accept values F^ = c2 < 0, for 
vacuum normalization reasons, and we are restricted to 
the upper half of the first quadrant. 

Four different regions show up. They are illustrated in 
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Fig. [S] and are characterized by the following modes 



A : 




^k^^, C^k'^^r^-', 












-9(1+3 




AB : 














- k-'^'+^-q-^", Cb ~ k^""?]-^- 


29 


(99) 


Bi : 
















29 


(100) 


B2 ■■ 




' k-'''+W~^\ C ^ k-'^'+W 


-29 










2q^ 


(101) 



The borders between the different regions are 

A~AB: cl = ^ 
AB-Bi: 4 



5 
3 

AB-B2: cl = 2- 
1 

9' 



2ci 



-B2 



(102) 

(103) 
(104) 
(105) 



The region where 'J is dominated by the spectrum of 
the first component is enlarged towards > 1. This 
include regions A and AB. In regions Bi and i?2, ^ gets 
the spectrum of the secondary component multiplied by 
fc^. In Si, perturbations are adiabatic and in B2 they 
are isocurvature. It is interesting to note how region Bi 
has been reduced by the cancellation discussed above. 
Only in region A do Cb and '^b take on the spectrum of 
the dominant component, while in all other regions they 
keep the original spectrum multiplied by k^. In region 
AB, therefore, we have coexistence of the two spectra in 
Ca and Cb, though C, is dominated by Ca- This is an effect 
of the other cancellation discussed before. 



VI. NUMERICAL SOLUTIONS 



too slow to be employed for an extensive investigation of 
the parameter space. As a good compromise, we use the 
default Mathematica method far from the bounce, where 
it is faster than Runge-Kutta and as reliable. Only in the 
interval |?7| < 10^^ we switch to a Runge-Kutta method 
which allows us to cross the bounce. After the bounce 
we return to the default Mathematica method until hori- 
zon re-entry. The combination of the two methods then 
gives a very efficient and reliable numerical integrator for 
perturbations in bouncing cosmologies. 

Since we are mainly interested in extracting the domi- 
nant mode at horizon re-entry, we have followed the evo- 
lution of two solutions with different k and we have eval- 
uated the spectrum of the generic variable / in the set 
{*,C,*b,Cb} as 



l0g/fci(??re) - l0g/fc,(?7r 

log fci - log k2 



(106) 



where the subscripts 1 and 2 refer to the two solutions 
with different values of fc, and rjre has been chosen of the 
order of 1/fc but small enough to avoid contamination 
from the first oscillation. Of course, using more solu- 
tions would have given a more accurate and reliable fit of 
the spectrum, but it turns out that the accuracy of the 
simulation allows this simple and fast choice. 

The other important prediction is on the time- 
dependence of the dominant mode, i.e. the exponent of 77 
(hereafter we shall call it time index) . This can be easily 
extracted as 



dlog/ 



dlogrj 



(107) 



v=nr 



We have tested two cosmological models: a bounce 
with two perfect fluids, and a bounce with a scalar field 
endowed with a tracking potential plus a perfect fluid. 
In both cases the secondary fluid has negative energy 
density in order to trigger the bounce. We shall present 
the results separately. 



A numerical test of the above predictions has been per- 
formed by integrating the set of differential equations 
(|43|I - H46(I using Mathematica. Initial conditions have 
been posed far from the bounce using the vacuum nor- 
malization prescription described in Sect. IIII Dl The de- 
fault Mathematica method for resolution of ordinary dif- 
ferential equations (which is a combination of Gear and 
Adams methods) cannot go through the bounce, because 
the step size falls to zero asymptotically as 77 — > 0. This 
is a consequence of the presence of 7i~^ factors in the dif- 
ferential equations, which requires a high precision can- 
cellation in the numerators. However, all variables and 
their derivatives have regular limits at the bounce, as can 
be explicitly checked on the so-obtained pre-bounce solu- 
tions. Thus a fixed step size method would overcome this 
problem. In fact, a standard Runge-Kutta method nicely 
follows the evolution of perturbations in any phases and 
in particular across the bounce. However, this method is 



A. Numerical testing of the two-perfect-fiuid model 

With this class of numerical simulations we can test 
the predictions of Sect. IV BI We have spanned the plane 
(c^,c^) starting from = 0.1 up to = 2 with 
running from -|- 0.1 up to — 2.5, with steps of 0.1 in 
both parameters. 

After the evaluation of the spectral and time index for 
all variables, just to have a global view of the parame- 
ter plane, we can draw contour plots for each index and 
compare these with the diagram of Fig. |S1 These are 
shown in Fig. |^ 

The orientation of the contours in the plots of Fig. 
d evidently shows the existence of two main regions of 
the parameter plane where the spectrum and the time 
dependence of the dominant mode in ^' and C at re- 
entry are different. The calculated border line between 
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FIG. 6: On the left panels we have the contour plots of the spectral indices and on the right panels the time indices in the 
case Ua = Ub = 0. (a) and (b) refer to (c) and (d) to f, (e) and (f) to ^i,, (g) and (h) to C^. 
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region AB and region B2 is superimposed in these fig- 
ures. Actually, the numerical border line seems to be 
slightly higher towards larger values of c^. However, this 
border line is very sensitive to the choice of the evalua- 
tion time rjre. Taking different times we can see that the 
numerical border line tends to the theoretical one in the 
limit rjre ^ l/k. In practice, one should keep in mind 
that close to the border line between regions AB and B2 
the two competing modes are comparable, giving rise to 
mixed initial conditions for the evolution inside the hori- 
zon. This also provokes glitches in the evaluation of the 
indices in these intermediate cases, which show up in our 
plots with a large density of contours on the upper side 
of the border line. 

In Fig. Et-h, the contours of spectral and time indices 
of and Cf) are represented together with the border line 
of regions A and AB. With the same caveat as before, 
also this theoretical border line is respected. 

Now, let us make a closer comparison between the the- 
oretical values of the indices and those found numerically. 
In Fig. Owe represent the results for a section of the plane 
(c^, cl) at fixed = 0.1 and varying cf. This line crosses 
three of the four theoretical regions, namely regions A, 
AB and Bi . The numerical values are plotted along with 
the lines representing the theoretical behaviour of the in- 
dices. We can see that the agreement is very satisfactory 
for all gauge-invariant variables. 

The spectra of \1/ and C are perfectly reproduced both 
in regions A, AB and in region Bi , where they are dom- 
inated by the spectrum of the secondary fluid. The dis- 
tinction between regions Bi and B2 is only in the time 
indices of 5* and (. In Figs. [7]p,d we have also plotted as 
a dashed line the theoretical behaviour that these vari- 
ables would have if region Bi did not exist and there were 
only region B2. We see that while the numerical values 
change but are still compatible with the expectation of 
region Bi , they cannot be compatible with the behaviour 
predicted by region i?2. 

The secondary fluid variables clearly show the transi- 
tion between regions A and AB. Only one point in "ifc 
still seems to follow the dominant mode of region A rather 
than that of region AB for cl = 0.7. As already stated, 
this depends on a conservative choice of Ty^e, which elimi- 
nates contaminations from the oscillations but can some- 
times select a not-yet-decayed mode in almost marginal 
cases like this. 

The existence of several modes that are of the same 
order at rjre is particularly evident for c^ ^ 1 and > 
1. In Fig. IHlwe have plotted the numerical values for 
the time index of 4* obtained for different choices of 77^6, 
along the line in the parameter plane with — 0.7. A 
conservative choice such as rjre — 0.4/(cafc) leads to very 
good agreement with the theoretical values far from the 
border line, but the points with of close to the border line 
between regions AB and B2 tend to follow the behaviour 
of region AB. However, pushing Ty^-e towards the first 
oscillation gives the necessary time to these points to feel 
the mode of region B2 and align to the others. However, 



the first oscillation globally alters all evaluations, which 
are displaced more and more from the theoretical super- 
horizon values. This discussion evidently shows that close 
to the border line the two competing modes contribute to 
determine the initial oscillations after horizon re-entry. 

B. Numerical testing of the scalar-field model 

In the second model of the bounce we have tested, the 
dominant fluid is replaced by a scalar field with an ex- 
ponential potential. As shown in Ref. , the exponent 
of the potential determines Ta in the range [0; 1]. At the 
same time, is fixed to 1. With this model, we can thus 
test the theoretical predictions for the case a 0. How- 
ever, the Fa > 1 part of the plane, which corresponds to 
region B2 of Fig. ^ remains unexplored numerically. One 
more limitation arises at small , where viable bouncing 
backgrounds cannot be obtained. A complete study of 
this class of backgrounds has been done only in the case 
Cfe = 1 E3> so that the reasons why the dynamics does 
not flow in the same way for small remains unknown 
to us. However, since we are mainly interested in per- 
turbations, we content ourselves with higher values of c^. 
This will be sufficient to distinguish regions Ai, A2 and 
Bi, as we shall see. 

It must also be said that in this case we have found 
it convenient to switch on a different set of differential 
equations close to the bounce. This set is the same as 
Eqs. (|43|I - H46(I but with a <-> &. In practice, to reduce 
numerical errors close to the bounce it is better to fol- 
low ^P, C, and (a- In this way, the whole numerical 
simulation is split in four phases: far pre-bounce, near 
pre-bounce, near post-bounce and far post-bounce. 

Figure |51 is the analogue of Fig. [7| for the bouncing 
cosmologies where the dominant component is a scalar 
field. However, here, we have managed to find good back- 
grounds only starting from = 0.5. In all the spectra 
(Figs. [7^,c,e,g) the transition from region A2 to region 
Bi at = 1, as predicted in Sect. IV AI is evident. Only 
one point testifies the existence of region Ai, where we 
see that spectrum and time dependence of and (b 
change. 

VII. CONCLUSIONS 

In this paper, we have addressed the problem of the 
evolution of perturbations in bouncing cosmologies by 
testing explicit regular models through the use of Ein- 
stein equations. At the moment, this is the only effec- 
tive and conclusive way to carry out a theoretical in- 
vestigation of the problem. In fact, in the absence of 
clean model-independent general arguments supporting 
any claim, the only answer can come from the study of 
regular bouncing models within the frame of General Rel- 
ativity, where the evolution of perturbations can be fol- 
lowed explicitly from the initial state to the final horizon 
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FIG. 7: On the left panels we have the spectral indices and on the right the time indices for our variables in the case a = 0, 
keeping = 0.1 and varying c^. (a) and (b) refer to ^, (c) and (d) to (e) and (f) to ^'6, (g) and (h) to Cd- 
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FIG. 8: Time index for * for cl = 0.7 and varying cf. The 
evaluation has been done at different rjre- Points along the 
solid line are taken at rjre ~ 0.4/(cafe). The dotted line is 
obtained with 77^6 ~ 1.5/(cafc), the dot-dashed line with rjre — 
2.5/{cak), and the dashed line with rjre = S.5/{cak). 

re-entry in the current expansion era. 

Within classical general relativity, a bounce from con- 
traction to expansion in spatially flat cosmologies can be 
realized only by violating the NEC. All high-energy cor- 
rections to General Relativity can thus be represented 
by a negative energy component in the Friedman equa- 
tions, which is important only in the bounce phase. With 
this philosophy, we have completed a full study of a 
very wide class of bouncing cosmologies, with two non- 
interacting components, allowing the primary to have an 
intrinsic isocurvature mode, in order to describe scalar 
fields and perfect fluids at the same time. We have 
shown that the possibly scale-invariant spectrum of ^' 
is carried by a homogeneous solution of the differential 
equations, which decays in the post-bounce phase. No 
other mode receives this spectrum, which is doomed to 
be only a transient characteristic of the Bardeen poten- 
tial. If we use vacuum-normalized initial conditions for 
both components, we have explicitly shown that this de- 
caying mode becomes subdominant w.r.t. some other 
mode before horizon re-entry. This general result ap- 
pears to contradict some claims made in the literature, 
in particular those of Refs. and ^3 ■ 

As already explained in our previous paper i in R^f. 
[Tl|. Eq. (26) provides an analytical solution /(r;) — 
77/(1 -|- 77^)^ for perturbations outside the horizon, which 
is valid both far from and through the bounce. However, 
the authors use an approximate form {fappiv) = ^/if) 
of this solution far from the bounce and then match fapp 
to / just before the bounce. Had they used / from the 
beginning throughout the bounce, they would have found 
no mixing. No surprise, therefore, that the mixing they 
find is at order 77""", which is just the next to leading order 
term in 7(77), which is missing in fappirj). Had they kept 
the r]~^ term in fappiv) they would have found mixing 
at the order 77^^. We thus conclude that the mixing 
they find is just an artifact of their matching procedure. 
We also note that a recent paper JJJ has criticized the 



splitting of ^ proposed by Ref. j^J, into two separate 
variables '5+ and which satisfy the same equation of 
motion as '5. 

In the case of Ref. the limit 7 ^ (fc ^ in 

our notation), used to evaluate the final spectrum, just 
selects the infrared modes at the bounce. However, the 
final spectrum should be evaluated at the horizon re- 
entry, after the decaying mode has gone away; it would 
correspond to 7 — > while keeping jz (our ^77) of order 
1. 

The most important lesson from this class of models 
is probably contained in the simplified example of Sect. 
IIVDI In fact, one may hypothesize any type of initial 
conditions for the perturbations of the secondary compo- 
nent, but one should then accept the basic facts, which 
are already included in the decoupled evolution of the 
perturbations of the dominant component. It should be 
concluded that in bouncing cosmologies with fast evolu- 
tion (g > 1/2) the decaying mode of the Bardeen poten- 
tial is always subdominant with respect to its constant 
mode, carrying a blue spectrum. In bouncing cosmolo- 
gies with slow evolution (g < 1/2) the two modes are of 
the same order at horizon re-entry. But then, even in 
the case when no other isocurvature mode coming from 
coupling to other sources dominates, the late-time phe- 
nomenology would be different from the observed one 

m 

There are at least three interesting directions in which 
this research may be extended. 

The first is the study of more complicated models 
for perturbations, going beyond the relations (Pn|) - (|^ . 
High-energy corrections may include higher derivatives 
and anisotropic stress, which are definitely outside the 
relations at the basis of our investigation. However, a 
first negative result already comes from Ref. [l^ , whose 
conclusions are similar to ours. 

The second direction is in allowing an interaction or a 
correlation between the two components. A more realis- 
tic model, also closer to the original ekpyrotic proposal, 
should include the decay of the dominant component to 
radiation at the bounce. The implementation of such 
models with a coherent treatment of classical perturba- 
tion does not seem trivial. 

Finally, the bounce integral may depend on an addi- 
tional scale of the theory, which is independent of the 
curvature scale of the bounce. It is difficult to imagine a 
physical component with these characteristics. However, 
the possibility remains open. 

In conclusion, all these possible extensions must in- 
voke some special physics at the bounce in order to up- 
set the basic result that the scale-invariant spectrum of 
the Bardeen potential is carried only by a decaying mode 
in the post-bounce. Up to now, no explicit construction 
of such a possibility has been presented in the literature 
and thus the claim that a scale- invariant spectrum can be 
carried by a constant Bardeen potential after a bounce 
remains an unproved conjecture. 
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FIG. 9: On the left panels we have the spectral indices and on the right the time indices for our variables in the case a ^ 0, 
keeping = 0.1 and varying . (a) and (b) refer to (c) and (d) to (e) and (f) to ^'a, (g) and (h) to ^b- 
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